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Abstract. Multiple orthogonal polynomials are a generalization of orthogonal polyno- 
mials in which the orthogonality is distributed among a number of orthogonality weights. 
They appear in random matrix theory in the form of special determinantal point pro- 
cesses that are called multiple orthogonal polynomial (MOP) ensembles. The correlation 
kernel in such an ensemble is expressed in terms of the solution of a Riemann-Hilbert 
problem, that is of size {r + 1) x (r + 1) in the case of r weights. 

A number of models give rise to a MOP ensemble, and we discuss recent results on 
models of non-intersecting Brownian motions, Hermitian random matrices with external 
source, and the two matrix model. A novel feature in the asymptotic analysis of the latter 
two models is a vector equilibrium problem for two or more measures, that describes 
the limiting mean eigenvalue density. The vector equilibrium problems involve both an 
external field and an upper constraint. 
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1. Introduction 

1.1. Random matrix theory. The Gaussian Unitary Ensemble (GUE) is 
the most prominent and most studied ensemble in random matrix theory. It is a 
probability measure on n x n Hermitian matrices for which the joint eigenvalue 
probability density function (p.d.f.) has the explicit form 

1 " 
Y- n i^k-x,rl[e~9^'^ (1) 

" l<j<k<n j=l 

where Z„ is an explicitly known constant. The density ([T|) can be analyzed with 
the help of Hermite polynomials. Due to this connection with classical orthogonal 
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polynomials many explicit calculations can be done, both for finite n and in the 
limit n — ^ oo, see |40| . In particular it leads to a description of the limiting 
behavior of eigenvalues on the global (macroscopic) scale as well as on the local 
(microscopic) scale. The global scale is given by the well-known Wigner semi-circle 
law 



p(x) = -L^4_a;2 -2<a;<2, 



(2) 



in the sense that for eigenvalues xi, . . . ,Xn taken from ([Ij, the empirical eigenvalue 
distribution ^ Sj=i '^(•^i) converges weakly to p{x) almost surely as n — )> oo. 
The local scale is characterized by the sine kernel 



S{x,y) 



sin tt{x — y) 
Tr(x - y) 



(3) 



in the bulk. This means that for any given x* G (—2, 2) and any fixed m e N, the 
TO-point correlation function (i.e., the marginal distribution) 



-Rm,n(^l, ■ . ■ , ^m) 



[n — m)\ 
has the scaling limit 



T n ( 

" T-<j<k<n 



dXm+l ■ ■ ■ dXn (4) 



lim \ 

n->oo [p[x*)n\ 



-Rn 



X + 



Xi 



p{x*)n 



. . , a; -I- 



p{x*)n 



= dei[S{x,,Xj\^,^^^^^. (5) 
At the edge points ±2 the sine kernel ([3]) is replaced by the Airy kernel 

Ai(:E)Ai'(y)-Ai'(a;)Ai(y) 



A{x,y) = 



x-y 



(6) 



and a scaling limit as in ([5]) (with scaling factor cn?^^ instead of p{x*)n) holds 
for X* — ±2. This result leads in particular to the statement about the largest 
eigenvalue 



lim Prob max Xk < 2 

n— >oo \ l<k<n 



,2/3 



det [I - A(t,oo)] 



(7) 



where A is the Airy kernel ([5]) and the determinant is the Fredholm determinant of 
the integral operator with Airy kernel acting on L^{t, oo). The limiting distribution 
([7]) is the famous Tracy- Widom distribution named after the authors of the seminal 
work |45] in which the right-hand side of ([7]) is expressed in terms of the Hastings- 
McLeod solution of the Painleve II equation. 
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These basic results of random matrix theory have been extended and gener- 
alized in numerous directions. Within the theory of random matrices, they have 
been generalized to ensembles with unitary, orthogonal and symplectic symmetry 
and to non- invariant ensembles (Wigner ensembles). The distribution functions of 
random matrix theory also appear in many other probabilistic models that have 
no apparent connection with random matrices (models of non-intersecting paths, 
tiling models, and stochastic growth models), see e.g. [8], [31] . 

Mehta's book [40] is the standard reference on random matrix theory. The book 
of Deift [22] has been very influential in introducing Riemann-Hilbert techniques 
into the study of random matrices. In recent years, a number of new monographs 
appeared [5], [B], [TB], [13], [H] that cover the various aspects of the theory of 
random matrices. 

1.2. Unitary ensembles and orthogonal polynomials. One direc- 
tion within random matrix theory is the study of ensembles of the form 

— e-"TrV(M)^^ (8) 

Zn 

defined on ri x n Hermitian matrices M, which reduces to the GUE in case V(x) = 
ia;^. The ensembles (|S|) have the property of unitary invariance and are called 
unitary ensembles. The eigenvalues have the p.d.f. 

j- n (-'c-x,)^ne-"^(^^) (9) 

" l<j<k<n j=l 

with a different normalizing constant Z„. [Throughout, we use Z„ to denote a 
normalizing constant, which may be different from one formula to the next.] 

Again explicit calculations can be done due to the connection with orthogonal 
polynomials [23] , [40] . For a given n, we consider the monic polynomial Pk^n of 
degree k that satisfies 



Pfc,„(x)x^e-"^(^)dx = hk,nSj,k, j=0,...,k. 
Then ((9)) is a determinant al point process [2], [44], with kernel 



Kn{x,y) = Ve-"^(^)Ve-"^fa)^ (10) 

fc=o 

which means that for every m G N the m-point correlation functions, defined as in 
(|4]), have the determinantal form 

det [K„{x,,Xj)]^^^^^^^^^^^^. 

As n — ?> cx), the limiting mean eigenvalue density 

p{x) — lim —Kn{x,x) 

n— 7-00 fi 
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is no longer Wigner's semi-circle law but instead it is the density p of the 
probability measure /i that minimizes the weighted logarithmic energy 



among all probability measures on JR. 

Local eigenvalue statistics, however, have a universal behavior as n — oo, that 
is described by the sine kernel ^ in the bulk. Thus for points x* with p{x*) > 
the limit ([5]) holds true. At edge points of the limiting spectrum the density p 
typically vanishes as a square root and then the universal Airy kernel ^ appears. 
For real analytic potentials V this was proved in [TT], [53] using Riemann-Hilbert 
methods. This was vastly extended to non-analytic potentials in recent works of 
Lubinsky [3 8) and Levin and Lubinsky |37] . among many others. 

1.3. This paper. In this paper we present an overview of the work (mainly of 
the author and co-workers) on multiple orthogonal polynomials and their relation 
to random matrix theory. Multiple orthogonal polynomials are a generalization of 
orthogonal polynomials that have their origins in approximation theory (Hermite- 
Pade approximation), see e.g. [3} l42]. 

They enter the theory of random matrices via a generalization of © which we 
call a multiple orthogonal polynomial (MOP) ensemble [34]. We present a number 
of models that give rise to a MOP ensemble, namely the model of non-intersecting 
Brownian motions, the random matrix model with external source and the two 
matrix model. 

The MOPs are described by a Riemann-Hilbert problem that may be used for 
asymptotic analysis as n — ^ cx) by extending the Deift-Zhou method of steepest 
descent |25| . The extensions are non-trivial and involve either an a priori knowledge 
of an underlying Riemann surface (the spectral curve) or the formulation of a 
relevant equilibrium problem from logarithmic potential theory [43], which asks 
for a generalization of the weighted energy functional ((TT]) . 

The latter approach has been succesfully applied to the random matrix model 
with external source and to the two matrix model, but only in very special cases, 
as will be discussed at the end of the paper. 



2. Multiple orthogonal polynomials 

2.1. MOP ensemble. We will describe here multiple orthogonal polynomi- 
als of type II, which we simply call multiple orthogonal polynomials. There is also 
a dual notion of type I multiple orthogonal polynomials. 

Suppose we have a finite number of weight functions wi, . . . ,Wr on M and 
a multi- index ft = (ni, . . . ,nr) G N^. Associated with these data is the monic 
polynomial Pa of degree \n\ — ni + ■■■+ rir so that 




(11) 




— oo 



Pfi{x)x^ Wk{x) dx = 0, for j = 0, . . . , 



rifc - 1, fc=l,...,r. 



(12) 
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The linear system of equations (|12l) may not be always uniquely solvable, but in 
many important cases it is. If Pfi uniquely exists then it is called the multiple 
orthogonal polynomial (MOP) associated with the weights Wi, . . . ,Wr and multi- 
index n. 

Existence and uniqueness does hold in the following situation. Assume that 



n (^'^ " 

l<j<fc<n 



(13) 



is a p.d.f. on R", where n = \n\ and the linear span of the functions /i, . . . , /« is 
the same as the linear span of the set of functions 

{x^Wkix) I j = 0, . . . - 1, fc = 1, . . . ,r}. 



So the assumption is that (|13p is non- negative for every choice of xi, . . . , x„ £ R", 
and that the normalization constant Z„ can be taken so that the integral (jl3p over 
R" is equal to one. Then the MOP satisfying exists and is given by 



Pftix) = E 



We call a p.d.f. on R" of the form ^ a MOP ensemble, see [31]. 

2.2. Correlation kernel and RH problem. The MOP ensemble (fT3l) 
is a determinantal point process [44j (more precisely a biorthogonal ensemble ^17; ) 
with a correlation kernel Kn that is constructed out of multiple orthogonal poly- 
nomials of type II and type I. It is conveniently described in terms of the solution 
of a Riemann-Hilbert (RH) problem. This RH problem for MOPs [47] is a gen- 
eralization of the RH problem for orthogonal polynomials due to Fokas, Its, and 
Kitaev W. 

The RH problem asks for an (r + 1) x (r -|- 1) matrix valued function Y so that 

• y : C \ R ^ ^("^+1)^ is analytic, 

• Y has limiting values on R, denoted by Y-^ and y_, where Y±{x) is 
the limit of Y{z) as z — > .t e R with ± Im z > 0, satisfying 

Wrix)\ 



Y+{x)^Y^ix) 



1 1 wi{x) 
1 



\0 

Y{z) = (/ + 0(l/z))diag(z" 







1 / 



(14) 



for X e 



-Ml 



z 



as z 



oo. 



If the MOP Pfi with weights wi, . . . ,Wr and multi-index ft = (rii, . . . , n^) exists 
then the RH problem (fT4| has a unique solution. If the MOPs with multi-indices 
n — Bj also exist, where ej is the jth unit vector of length r, then the first column 
of Y consists of 

Yi,i(z) =P«(z), 



Yj + l,l{z) = Cj^nPn~ej (z), j = 1, . . . , r 



(15) 
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where Cj^s is the constant 



-27ri 



Ppi-Sj (a;)a;"^ Wj {x)dx 



^0. 



The other cohimns of Y contain Cauchy transforms 



1 



Yjs{x)wk {x) 



dx. 



, ,r + 1, = 1, . 



It is a remarkable fact that the correlation kernel of the MOP ensemble (|13p is 
expressed as follows in terms of the solution of the RH problem, see j20| . 



1 



2TTi{x — y) 



(0 w,{y) ■■■ Wr{y))Y^\y)Y+{x) 







x,y e 



(16) 

The inverse matrix F ^ contains MOPs of type I, and the formula is essentially 
the Christoffel-Darboux formula for multiple orthogonal polynomials. 

Besides giving a concise formula for the correlation kernel, the expression (jl6p 
for the kernel gives also a possible way to do asymptotic analysis in view of the 
Deift-Zhou method of steepest descent for RH problems. 



3. Non intersecting path ensembles 

A rich source of examples of determinantal point processes is provided by non- 
intersecting path ensembles. In special cases these reduce to MOP ensembles. 

3.1. Non-intersecting Brownian motion. Consider a one-dimensional 
strong Markov process with transition probability densities Ptix, y) for t > 0. Sup- 
pose n independent copies are given with respective starting values oi < 02 < 
• • • < a„ at time t = Q and prescribed ending values &i < &2 < • • • < at time 
t — T > Q that are conditioned not to intersect in the full time interval G <t <T. 
Then by an application of a theorem of Karlin and McGregor [32] , the positions 
of the paths at an intermediate time t G (0, T) have the joint p.d.f. 

^ det [pt{aj,Xk)]i^j^k<n ■ det [PT-t(a;fc, 60]i<fe,i<„ • (17) 

In a discrete combinatorial setting the result of Karlin and McGregor is known as 
the Lindstrom-Gessel-Viennot theorem. 

The density function (|17|) is a biorthogonal ensemble, which in very special 
cases reduces to the form ([T3|) of a MOP ensemble. 
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Figure 1. Non-intersecting Brownian bridges starting and ending at 0. At any inter- 
mediate time t £ (0, 1) the positions of the paths have the same distribution as the 
(appropriately rescaled) eigenvalues of an n x n GUE matrix. 



An example is the case of Brownian motion (actually Brownian bridges) with 
the transition probability density 

, , 1 (^-y)^ 

Pt{x,y) = ^— e 2t , t>Q. 



^ n 



In the confluent limit where all Uj — > and all 6; — > the p.d.f. (|17l) turns into 
1 

^" l<j<k<n 3 = 1 

with a different constant Z„. This is up to trivial scaling the same as the p.d.f. (H]) 
for the eigenvalues of an n x n GUE matrix. 

If however, we let all 6; — and choose only r different starting values, denoted 
by fli, . . . , Or, and Uj paths start at aj, then (|17p turns into a MOP ensemble with 
weights 

u;,(2;) = e~^^(^"'+^", j = (18) 

and multi-index (ni, . . . , rir). This is a multiple Hermite ensemble, since the asso- 
ciated MOPs are multiple Hermite polynomials [5] 

3.2. Non-intersecting squared Bessel paths. The squared Bessel 
process is another one-dimensional Markov process which gives rise to a MOP 
ensemble. The squared Bessel process is a Markov process on [0,oo), depending 
on a parameter a > — 1, with transition probability density 

P*(.,y)==^©"'%--^-^)/J^), .,.>o. 
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Figure 2. Non-intersecting Brownian bridges starting at two different values and ending 
at 0. At any time t G (0, 1) the positions of the paths have the same distribution as the 
eigenvalues of an n x n GUE matrix with external source. The distribution is a multiple 
Hermite ensemble with two Gaussian weights (|18p . 

where la is the modified Bessel function of first kind of order a. In the limit where 
all Qj ^ a > and bj — ^ the p.d.f. ([T7| for the positions of the paths at time 
t e (0, T) is a MOP ensemble with two weights 



and multi-index (rii,ri2) where ni — \n/2~\ and n2 — ln/2\, see |35j . In the limit 
a — > this further reduces to an orthogonal polynomial ensemble for a Laguerre 
weight. 



4. Random matrix models 

The random matrix model with external source, and the two matrix model also 
give rise to MOP ensembles. 

4.1. Random matrices with external source. The Hermitian matrix 
model with external source is the probability measure 




J_g-n Tr(y (M) -AM) 



(19) 



Multiple orthogonal polynomials 



9 



on n X n Hermitian matrices, where the external source A is a given Hermitian 
n X n matrix. This is a modification of the usual Hermitian matrix model, in which 
the unitary invariance is broken (TS], [15] . 

Due to the Harish-Chandra/Itzykson-Zuber integral formula 30 , it is possible 
to integrate out the eigenvectors explicitly. In case the eigenvalues ai, . . . , a„ of A 
are all distinct, we obtain the explicit p.d.f. 

1 " 

— det[e--^],<,_^.<„. n {xk-x,)-l[e-''^^^^ 

" l<j<k<n j=l 

for the eigenvalues of M. In case that oi, . . . , are the distinct eigenvalues of A, 
with respective multiplicities ni, . . . , n^, then the eigenvalues of M are distributed 
as a MOP ensemble ([T5)) with weights 

= e-"(^(^)-'^^-^\ j = l,...,r (20) 

and multi- index (ni, ...,71^), see |13] 

For the case V{x) = \x'^ the external source model is equivalent to the 
model of non-intersecting Brownian motions with several starting points and one 
ending point, cf. (fT8|) . 

4.2. Two matrix model. The Hermitian two matrix model 

_}_^-mr(V(M^) + W(M.)-rM,M,) ^^^^^2 (21) 
Zn 

is a probability measure defined on couples (Mi, M2) oi nx n Hermitian matrices. 
Here V and W are two potentials (typically polynomials) and r ^ is a coupling 
constant. The model is of great interest in 2d quantum gravity [21], [30], [33], as 
it allows for a large class of critical phenomena. 

The eigenvalues of the matrices Mi and M2 are fully described by biorthogonal 
polynomials. These are two sequences {Pk,n)k and {Qj^n)j of monic polynomials, 
degPfc,„ = fc, degQjVn = j, such that 

/oo /'OO 
/ Pfc,„(a:)g,,„(y)e-"(^(-)+'^(^)-^-^)dxdy = 0, ifj^fc, (22) 

see e.g. [9], [27], [40], [41]. 

If is a polynomial then the biorthogonality conditions ([22]) can be seen as 
multiple orthogonal polynomial conditions with respect to r = deg W — 1 weights 

/oo 
yOe-"(^(y^-^='y^dy, j = 0,...,r-l, (23) 
-00 

see [55|. Furthermore, the eigenvalues of Mi are a MOP ensemble with the 
weights (j23p and multi-index ft = (no, . . . , Wr-i) with = [n/r] for j = 0, . . . , q—1 
and rij — [n/rj for j = q, . . . , r — 1 if n = pr + q with p and < g < r non- negative 

integers, see [26j for the case where W{y) = 
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Figure 3. Non-intersecting Brownian bridges starting at two different values and ending 
at 0. As n — >■ oo, the paths fiU out a heart-shaped domain. Critical behavior at the cusp 
point is desribed by the Pearcey kernel (|24|l . 

5. Large n behavior and critical phenomena 

We discuss the large n behavior in the above described models. 

5.1. Non-intersecting Brownian motion. In order to have interesting 
limit behavior as n — > oo in the non-intersecting Brownian motion model we scale 
the time variables T !—> l/n, t i-> t/n, so that < t < 1. In the case of one 
starting value and one ending value, see Figure [U the paths will fill out an ellipse 
as n — >■ oo. 

In the situation of Figure [5] the paths fill out a heart-shaped region as n — ^ oo, 
as shown in Figure |3l New critical behavior appears at the cusp point where the 
two groups of paths come together and merge into one. 

Around the critical time the correlation kernels have a double scaling limit, 
which is given by the one-parameter family of Pearcey kernels 



where p and q are solutions of the Pearcey equations p"'(x) = xp(x) — hp'{x) and 
q"'{y) = yq{y) + bq'iv)- This kernel was first identified by Brezin and Hikami |18) 
who also gave the double integral representation 



V{x,y-h) 



P{x)q"{y) ~p'{x)q'{y) +p"{x)q{y) - hp{x)q{y) 



b e R, (24) 



x-y 




(25) 



where the contour C consists of the rays from zbooe*''/^ to and the rays from 
to iooe""/''. 
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Consideration of multiple times near the critical time leads to an extended 
Pearcey kernel and the Pearcey process given by Tracy and Widom |46l . 

As already noted above, the model of non-intersecting Brownian motion with 
two starting points and one ending point is related to the Gaussian random matrix 
model with external source 

1 g-nTr(iM^-AAf)^^^^ (26) 

with external source 

A — diag(a, . . . , a, —a, . . . , —a). (27) 

n/2 times n/2timcs 

In this setting the critical a- value is acrit — 1 and the Pearcey kernel ^1% arises 
as n -> oo with a = 1 + 1^1= ■ In. [15] this was studied with the use of the 



Riemann-Hilbert problem (jl4p for multiple Hermite polynomials with two weights 
g-n(ix ±ax)_ rpj^g asymptotic analysis as n — 00 was done with an extension of the 
Deift-Zhou method of steepest descent [IS] to the case of a 3 x 3 matrix valued RH 
problem. See also [M] and [4] for a steepest descent analysis of the RH problem 
in the non-critical regimes a > 1 and < a < 1, respectively. 

Another interesting asymptotic regime is the model of non-intersecting Brow- 
nian motion with outliers. In this model a rational modification of the Airy kernel 
appears that was first described in [7j in the context of complex sample covariance 
matrices, see also [Tj. 

5.2. Random matrices with external source. If V is quadratic in 
the random matrix model with external source ()19|) then this model can be mapped 
to the model of non-intersecting Brownian motions. Progress on this model beyond 
the quadratic case is due to McLaughlin [3£| who found the spectral curve for the 
quartic potential V(x) = \x'^ and for a sufficiently large (again A is as in ([27]) ). 

A method based on a vector equilibrium problem was introduced recently by 
Bleher, Delvaux and Kuijlaars [10 . The vector equilibrium problem extends the 
equilibrium problem for the weighted energy (jlip that is important for the unitary 
ensembles and which is crucial in the steepest descent analysis of the RH problem 
for orthogonal polynomials [24 . 

In jlOj it is assumed that V is an even polynomial, and that A is again given 
as in ()27p . The vector equilibrium problem involves two measures /ii and [12^ and 
it asks to minimize the energy functional 

log- — ^- — -dni{x)dfii{y) + 1 1 log- — - — Tdfi2{x)dfi2{y) 



\x-y\ J J \x~y\ 



log -dfii{x)d^2iy) + / {V{x) - a|a;|) dfii{x) (28) 

\x-y\ J 

where fii is on R with J d^i = 1, ^2 is on iM (the imaginary axis) with J diJL2 = 1/2, 
and in addition /i2 < cr, where a is the measure on iM. with constant density 

\dz\ TT 
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There is a unique minimizer, and the density pi of the measure /^i is the hmiting 
mean eigenvalue density 

pi{x) = hm -Kn{x,x) 

where Kn is the correlation kernel of the MOP ensemble with weights e~"(^(^)^°^). 
The RH problem (|14l) is analyzed in the large n limit with the Deift/Zhou steepest 
descent method in which the minimizcrs from the vector equilibrium problem play 
a crucial role. 

The upper constraint /i2 < f is not active for large enough a and in that case 
the support of /ii has a gap around 0. For smaller values of a the constraint a is 
active along an interval [— ?c, ic], c > 0, on the imaginary axis. Critical phenomena 
take place when either the constraint becomes active, or the gap around closes, 
or both. If one of these two phenomena happens, then this generically will be a 
phase transition of the Painleve II type that was described in the unitary matrix 
model in [12 and jl9j . If the two phenomena happen simultaneously then this 
is expected to be phase transition of the Pearcey type which, if true, would be a 
confirmation of the universality of the Pearcey kernels ([M)) at the closing of a gap 

m- 

Both kinds of transitions are valid in the external source model with even 
quartic potential V{x) — jx'^—^x'^, see [TU]. For the particular value t = -\/3, there 
is a passage from the Painleve II transition (for t > -\/3) to the Pearcey transition 
(for t < VS). The description of the phase transition for t = a/S remains open. 

5.3. Two matrix model. In [26] Duits and Kuijlaars applied the steepest 
descent analysis to the RH problem for the two matrix model (|21|) with quartic 
potential 

W(y) = y (30) 

and for V an even polynomial. The corresponding MOP ensemble has three weights 
of the form ([^5]) and the RH problem (|14l) is of size 4x4. Again a vector equilibrium 
problem plays a crucial role. 

The vector equilibrium problem in [26| involves three measures fJ-i^ fJ-2 and (j.^. 
It asks to minimize the energy functional 

3 



iff ^°8j;;;^^dfi,{x)dfij{y) 



- E// ^ogj-^^dp,{x)dp,+^{y) + J {V{x) - ^\Tx\^/')dp^{x) (31) 

among measures /ii on R with J dfii — 1, fi2 on iM. with J dii2 =2/3 and fi^ on 
K. with J dfi:i ==1/3. In addition fJ.2 l£ <y where cr is a given measure on iM. with 
density 
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There is a unique minimizer and the density pi of the first measure /xi is equal 
to the limiting mean eigenvalue density of the matrix Mi in the two matrix model. 
In addition, the usual scaling limits (sine kernel and Airy kernel) are valid in the 
local eigenvalue regime, see [26j . However there is no new critical behavior in the 
two matrix model with W is given by (pO| . 

New multicritical behavior is predicted in |21| for more general potentials. For 
the more general quartic potential W{y) = jy* — |?/^ an approach based on a 
modification of the vector equilibrium problem (j3ip is under current investigation. 
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